DischargeRoute

public subroutine DischargeRoute(dt, time, flowdir, runoff, dtrk, riverToGroundwater, groundwaterToRiver, storage)

route discharge in surface network

Arguments

Type	Intent		Name
integer(kind=short),	intent(in)	::	dt	time step [s]
type(DateTime),	intent(in)	::	time
type(grid_integer),	intent(in)	::	flowdir	flow direction
type(grid_real),	intent(in)	::	runoff	net rainfall [m/s]
integer(kind=short),	intent(in)	::	dtrk	time step for reservoirs
type(grid_real),	intent(in)	::	riverToGroundwater	discharge from river to groundwater (m3/s)
type(grid_real),	intent(in)	::	groundwaterToRiver	discharge from groundwater to river (m3/s)
type(grid_real),	intent(inout)	::	storage	volume stored on channel cell

Variables

Type	Visibility	Attributes		Name
real(kind=float),	public		::	Qdiverted	discharge diverted from a diversion channel (m3/s)
real(kind=float),	public		::	Qdownstream	discharge downstream a diversion channel (m3/s)
real(kind=float),	public		::	Qlateral	lateral discharge (m3/s)
real(kind=float),	public		::	QlateralChannel	lateral discharge in a diversion channel (m3/s)
real(kind=float),	public		::	area
type(Diversion),	public,	POINTER	::	d	pointer to current diversion
real(kind=float),	public		::	ddx
type(Diversion),	public,	POINTER	::	dnest	pointer to current diversion
integer(kind=short),	public		::	iin
integer(kind=short),	public		::	is
integer(kind=short),	public		::	jin
integer(kind=short),	public		::	js
integer(kind=short),	public		::	k
type(Reservoir),	public,	POINTER	::	p	pointer to current reservoir
real(kind=float),	public		::	runoff_ij	runoff of single cell
real(kind=float),	public		::	tWidth
real(kind=float),	public		::	wDepth

Source Code

SUBROUTINE DischargeRoute  &
!
(dt, time, flowdir, runoff, dtrk,  riverToGroundwater, &
    groundwaterToRiver, storage )

IMPLICIT NONE

!Arguments with intent(in):
INTEGER(KIND = short), INTENT(IN)    :: dt !!time step [s]
TYPE(DateTime), INTENT(IN)           :: time
TYPE(grid_integer), INTENT(IN)       :: flowdir !!flow direction
TYPE(grid_real), INTENT(IN)          :: runoff !! net rainfall [m/s]
INTEGER(KIND = short), INTENT(IN)    :: dtrk !!time step for reservoirs 
TYPE(grid_real), INTENT(IN)          :: riverToGroundwater  !!discharge 
                                            !!from river to groundwater (m3/s)
TYPE(grid_real), INTENT(IN)          :: groundwaterToRiver !!discharge 
                                             !!from groundwater to river (m3/s)
!Arguments with intent inout
TYPE(grid_real), INTENT(INOUT) :: storage !!volume stored on channel cell 

!local declarations:
INTEGER(KIND = short)    :: k, iin, jin, is, js
REAL (KIND = float)      :: ddx
REAL (KIND = float)      :: Qlateral !!lateral discharge (m3/s)
REAL (KIND = float)      :: QlateralChannel !!lateral discharge in a diversion channel (m3/s)
REAL (KIND = float)      :: Qdownstream !! discharge downstream a diversion channel (m3/s)
REAL (KIND = float)      :: Qdiverted !! discharge diverted from a diversion channel (m3/s)
REAL (KIND = float)      :: runoff_ij !!runoff of single cell 
REAL (KIND = float)      :: wDepth
REAL (KIND = float)      :: tWidth
REAL (KIND = float)      :: area
TYPE (Reservoir), POINTER :: p !!pointer to current reservoir
TYPE (Diversion), POINTER :: d, dnest !!pointer to current diversion
!-------------------------------end of declaration----------------------------- 

!reset Qin 
Qin  = 0.

!loop troughout hydro network
DO k = 1, streamNetwork % nreach
  iin = streamNetwork % branch (k) % i0
  jin = streamNetwork % branch (k) % j0
  
  !follow the branch
  DO WHILE (  .NOT.( jin == streamNetwork % branch (k) % j1 .AND. &
                    iin == streamNetwork % branch (k) % i1  )   )
    !set runoff_ij
    runoff_ij = runoff % mat(iin,jin)
    
     
    !find downstream cell
    CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), &
                        is,js, ddx, flowdir)
    
    ! looking for reservoir
    IF ( dtrk > 0) THEN
        IF ( dams % mat (iin,jin) /= dams % nodata ) THEN
            
            !set current reservoir
            p => GetReservoirById (list = pools, id = dams % mat (iin,jin) )
            
            IF (time == p % tReadNewStage) THEN
                !read new stage value from file
                CALL ReadData (network = p % network, &
                            fileunit = p % unit, time = time)
                
                !update target level
                p % stageTarget = p % network % obs (1) % value
                
                 p % tReadNewStage = p % network % time
                
                !update time of next reading of target stage
                p % tReadNewStage = p % tReadNewStage + &
                                    p % network % timeIncrement
            
            END IF
            
            
            IF (p % dischargeDownstream .AND. &
                time == p % tReadNewDischargeDownstream ) THEN
                !read new downstream discharge value from file
                CALL ReadData (network = p % networkDischargeDownstream, &
                            fileunit = p % unitDischargeDownstream, time = time)
                
                !update time of next reading of downstream discharge
                p % tReadNewDischargeDownstream = &
                         p % tReadNewDischargeDownstream + &
                         p % networkDischargeDownstream % timeIncrement
            END IF
                
            IF (p % dischargeDiverted .AND. &
                time == p % tReadNewDischargeDiverted ) THEN
                !read new diverted discharge value from file
                CALL ReadData (network = p % networkDischargeDiverted, &
                            fileunit = p % unitDischargeDiverted, time = time)
                
                !update time of next reading of diverted discharge
                p % tReadNewDischargeDiverted = &
                         p % tReadNewDischargeDiverted + &
                         p % networkDischargeDiverted % timeIncrement
            END IF    
           
            !reservoir routing
            Qin % mat(iin,jin) = Qin % mat(iin,jin) + &
                                    runoff_ij * &
                                    CellArea(runoff,iin,jin) 
            CALL LevelPool (time, dt, dtrk, Pin % mat(iin,jin), &
                            Qin % mat(iin,jin), p)
        
            !check and simulate diversion from reservoir
            IF ( p % bypassIsPresent ) THEN
                !divert flow from river
                d => p % bypass
        
                CALL DivertFlow (time, d, p % Qout, p % Qout )
                
                !overwrite diversion inlet discharge when observation is available
                IF ( p % dischargeDiverted ) THEN
                   IF ( p % networkDischargeDiverted % obs (1) % value /= &
                        p % networkDischargeDiverted % nodata ) THEN
                      d % QinChannel = p % networkDischargeDiverted % obs (1) % value 
                   END IF
                END IF
        
               !route discharge into channel
               Qlateral = 0.
               CALL MuskingumCungeTodini ( dt, &
                                d % channelLenght, &
                                d % QinChannel, &
                                d % PinChannel, &
                                d % PoutChannel, &
                                Qlateral, Qlateral, &
                                d % channelWidth, &
                                d % channelBankSlope, &
                                d % channelSlope, &
                                d % channelManning, &
                                d % QoutChannel, &
                                wDepth, tWidth )
           
               !copy data of current step for next step
               d % PinChannel  = d % QinChannel
               d % PoutChannel = d % QoutChannel
    
            END IF
           
            !set Qout
            SELECT CASE (  p % typ )
            CASE ( 'on' )
                 IF ( p % dischargeDownstream ) THEN
                   IF ( p % networkDischargeDownstream % obs (1) % value /= &
                        p % networkDischargeDownstream % nodata ) THEN
                      !overwrite reservoir downstream discharge
                      p % Qout = p % networkDischargeDownstream % obs (1) % value 
                   END IF
                END IF
                Qout % mat(iin,jin) = p % Qout
            CASE ( 'off' )  !off-stream reservoir  DA RIVEDERE!!
                !compute off-stream pool outflow
                CALL TableGetValue ( valueIn =  p % stage, tab = p % geometry, keyIn = 'h', &
                                    keyOut ='Qout', match = 'linear', &
                                    valueOut = p % Qout_off, bound = 'extendlinear' )
                p % Pout_off = p % Qout_off
                Qout % mat(iin,jin) = p % Qout
                !assign outflow to pool out cell
                !IF (p % rout /= MISSING_DEF_INT .AND. p % cout /= MISSING_DEF_INT ) THEN
                !Qin % mat(p % rout,p % cout) = Qin % mat(p % rout,p % cout) + &
                !    p % Qout_off
                !END IF
             
            END SELECT
            Pin % mat(iin,jin) = Qin % mat(iin,jin)
            Pout % mat(iin,jin) = Qout % mat(iin,jin) 
            jin = js
            iin = is
            CYCLE
        END IF
    END IF
        
    !Muskingum-Cunge-Todini
    
    area = CellArea(runoff,iin,jin)
        
    !set Qlateral
    Qlateral = runoff_ij * area
        
    !remove irrigation
    IF ( ALLOCATED (Qirrigation % mat) ) THEN
        Qlateral = Qlateral - Qirrigation % mat (iin, jin)
    END IF
        
    !include river-groundwater exchange
    IF ( ALLOCATED ( riverToGroundwater % mat ) ) THEN
        IF ( riverToGroundwater % mat (iin,jin) /= &
                riverToGroundwater % nodata ) THEN
           
            Qlateral = Qlateral + groundwaterToRiver % mat (iin,jin) - &
                                    riverToGroundwater % mat (iin,jin) 
            
        END IF
    END IF
    
    !add excess of water retained in snow
    IF ( ALLOCATED ( excessWaterSnowFlux % mat ) ) THEN
        Qlateral = Qlateral + excessWaterSnowFlux % mat (iin, jin) * area
    END IF
        
        
    !add outflow from by-pass channel or off-stream basin
    p => pools
    DO WHILE (ASSOCIATED(p)) !loop trough all reservoirs
        IF ( p % typ == 'off' ) THEN
            IF ( iin == p % rout .AND. jin == p % cout ) THEN
                    Qlateral = Qlateral + p % Pout_off
            END IF
        END IF
            
        IF (p % bypassIsPresent) THEN
            IF ( iin == p % bypass % rout .AND. jin == p % bypass % cout ) THEN
                    Qlateral = Qlateral + p % bypass % PoutChannel
            END IF
        END IF
            
        p => p % next
    END DO
    
    !remove diverted discharge from diversion channel
    IF (dtDiversion > 0) THEN
        IF ( divChannels % mat (iin,jin) /= divChannels % nodata ) THEN
            !divert flow from river
            d => GetDiversionById (list = diversionChannels, &
                              id = divChannels % mat (iin,jin) )
        
            CALL DivertFlow (time, d, Qin % mat(iin,jin), Qdownstream )
            
            Qdiverted = Qin % mat (iin, jin) - Qdownstream
           
            !update Qlateral
            Qlateral = Qlateral -  Qdiverted
            
            !route discharge into channel
            QlateralChannel = 0.
            CALL MuskingumCungeTodini ( dtDiversion, &
                                d % channelLenght, &
                                d % QinChannel, &
                                d % PinChannel, &
                                d % PoutChannel, &
                                QlateralChannel, QlateralChannel, &
                                d % channelWidth, &
                                d % channelBankSlope, &
                                d % channelSlope, &
                                d % channelManning, &
                                d % QoutChannel, &
                                wDepth, tWidth )
    
           !copy data of current step for next step
            d % PinChannel  = d % QinChannel
            d % PoutChannel = d % QoutChannel
        END IF
    END IF
    
    
        
    !add outflow from diversion channel
    d => diversionChannels
        
    DO WHILE (ASSOCIATED(d)) !loop trough all diversion channels
            
        IF ( iin == d % rout .AND. jin == d % cout ) THEN
            ! As the flood routing is solved from upstream to downstream
            ! PoutChannel contains the current or the previous time step discharge according to the 
            ! Horton order of the outlet section respect to the horton order of onlet section 
            Qlateral = Qlateral + d % PoutChannel
        END IF
            
        d => d % next
    END DO
        
    CALL MuskingumCungeTodini ( dt, ddx, Qin % mat(iin,jin), &
                                Pin % mat(iin,jin), &
                                Pout % mat(iin,jin), &
                                Qlateral, Plat % mat(iin,jin), &
                                sectionWidth % mat (iin,jin), &
                                bankSlope % mat (iin,jin), &
                                streamNetwork % branch (k) % slope, &
                                manning % mat (iin,jin), &
                                Qout % mat(iin,jin), &
                                waterDepth % mat(iin,jin), &
                                topWidth % mat(iin,jin) )
    
    storage % mat (iin,jin) = storage % mat (iin,jin) + &
                ( Qin % mat(iin,jin) + Qlateral - Qout % mat(iin,jin) ) * &
                dt / CellArea(Qin,iin,jin)
           
    IF ( storage % mat (iin,jin) < 0. ) THEN
        !storage % mat (iin,jin) = 0.
        !Qin % mat(iin,jin) = 0.
        !Qout % mat(iin,jin) = 0.
    END IF
              
    IF (isnan (Qout % mat(iin,jin)) ) THEN
        Qout % mat(iin,jin) = Pout % mat(iin,jin)
    END IF
      
    !computed Qout becomes Qin for the downstream section. 
    !Take account of junctions, sum to existing discharge
    Qin % mat(is,js) = Qin % mat(is,js) + Qout % mat(iin,jin)
       
    !store previous values for next time step
    Pin % mat(iin,jin) = Qin % mat(iin,jin)
    Pout % mat(iin,jin) = Qout % mat(iin,jin)
    Plat % mat(iin,jin) = Qlateral
    
    !store cell indexes and select downstream cell for next loop
    jin = js
    iin = is
    
    !check next cell is outlet cell
    CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), &
                        is,js, ddx, flowdir)
    
    IF ( CheckOutlet (iin,jin,is,js,flowdir) ) THEN
        !next cell will not be computed, set approximate value of Qout
        Qout % mat (iin, jin) =  Qin % mat (iin, jin)  + &
            runoff % mat(iin,jin) *  CellArea(runoff,iin,jin) 
    END IF
    
    
    
    
  END DO
END DO

!write results
IF ( time == timePointExport ) THEN
   CALL DischargePointExport ( time )
   timePointExport = timePointExport + sitesDischarge % timeIncrement
END IF

IF ( time == timePoolsExport ) THEN
   CALL OutReservoirs ( pools, time, Qin, Qout )
   timePoolsExport = timePoolsExport + dtOutReservoirs
END IF

IF ( time == timeDiversionsExport ) THEN
   CALL OutDiversions ( diversionChannels, time, Qin, Qout )
   timeDiversionsExport = timeDiversionsExport + dtOutDiversion
END IF


RETURN
END SUBROUTINE DischargeRoute

DischargeRoute Subroutine

Contents

Variables

Source Code

public subroutine DischargeRoute(dt, time, flowdir, runoff, dtrk, riverToGroundwater, groundwaterToRiver, storage)

Arguments

Variables

Source Code